home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / hips / sources / tools / mean_scale.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-08-06  |  20.7 KB  |  806 lines
  1. /*
  2. % MEAN_SCALE.C - scale the image to entire mean of background mean or
  3. %        weighted forground mean value or certain given value.
  4. %
  5. %    Copyright (c)    Jin Guojun
  6. %
  7. % compile:
  8. %    cc -O -o mean_scale mean_scale.c -lscs3 -lccs -lhips -lrle -ltiff -lm
  9. %
  10. % AUTHOR:    Jin Guojun - LBL    1/11/91                */
  11.  
  12. char    serbuf[256], usage[]="options\n\
  13. -a    adjust float input automatically.\n\
  14. -A x y [w [h]]    use given Area (center at x..y with size w,h [16x16]) as\n\
  15.     mean-scale reference.    \n\
  16. -B[#]    mean_scale to weight Balanced brightness.\n\
  17. -b[#]    make all background have same value. If further value given    \n\
  18.         immediately, then make background to that value.    \n\
  19. -C \"start_postion square_size\"    \n\
  20.     choose 4 corners, mean them as background value.    \n\
  21. -c    mean_scale to ceiling. default is mean_scale to floor.    \n\
  22. -e[#]    enhanced mode. # is enhance factor for trying on the worst case.\n\
  23.     No necessary to use it at regular time.    \n\
  24. -E #    End frames special. # < 0, 0..#; # > 0, n-#..n    \n\
  25. -F #    overflow value handle Factor. Default=1. If F=1.5, then when    \n\
  26.     overflow happened, that pixel value = max / F. It is used for    \n\
  27.     avoiding too bright in a certain case.        \n\
  28. -f[#]    foreground threshold. A value won't be changed if the value is    \n\
  29.     lower then this threshold.    \n\
  30. -h[#]    height percentage from the top of a hill. It'd better work with    \n\
  31.     -S option. It is defaulted to 50%.    \n\
  32. -l #    \n\
  33. -m #    \n\
  34.     thresholds -- see message at begin of running.    \n\
  35.     The -m is main threshold.    \n\
  36.     -l means that any value lower than -l value # are treated as 0.    \n\
  37.     It is used to build narrow hills.\n\
  38. -n # [#s]    num frames, start frame    \n\
  39. -S[#]    smoothing search. An enhanced method for precisely process.    \n\
  40.     The good range is 1-3. Default=2 (also for enhanced mode    \n\
  41.     without smooth option or smooth factor greater than 3.        \n\
  42. -s #    step (sample width). Default=3. This is a key parameter in    \n\
  43.         MEAN_SCALE processing. The width 2 or 4, sometimes,    \n\
  44.         works better than 3 for BYTE format image. So if result    \n\
  45.         is not desired, try -s 2 or 4 with/without -e options    \n\
  46.         before trying other options.    \n\
  47. -T    Testing for single frame.        \n\
  48. -t    count top values, otherwise omit.    \n\
  49. -V #    valley position. The first valley in main chart is used if -v    \n\
  50.     without -V.\n\
  51. -v[#]    dig valley with width #. Default=samples width.    \n\
  52. -W #    min hill foot width. default=2.5 sample width.    \n\
  53. -w #    min weight of choosing hill.\n\
  54. -z[#]    evaluate zeros. BUT if # given, then not evaluate all value below #\n\
  55. -G    send all information to a file which name is input file name    \n\
  56.     plusing all options and suffix \".rpt\".\n\
  57. -M    print all message on screen.        \n\
  58. -o outfile_name                    \n\
  59.     thie is used for PC kind computer which requires binary output    \n\
  60.     file mode.    \n\
  61. [<] in_filename [> out_file]\n";
  62.  
  63. #include "header.def"
  64. #include "imagedef.h"
  65. #include <math.h>
  66.  
  67. U_IMAGE    uimg;
  68.  
  69. #define    ibuf    uimg.src
  70. #define    obuf    uimg.dest
  71. #define    cln    uimg.width
  72. #define    row    uimg.height
  73. #define    pxl_bytes    uimg.pxl_in
  74.  
  75. #ifdef    IBMPC
  76. #define    MinReset    32767
  77. #else
  78. #define    MinReset    131072
  79. #endif
  80. #ifndef    SampleWidth
  81. #define    SampleWidth    3
  82. #endif
  83. #ifndef    MOSTHILLS
  84. #define    MOSTHILLS    16
  85. #endif
  86. #ifndef    THRESHOLD
  87. #define    THRESHOLD    32
  88. #endif
  89. #ifndef    VISIBLE_LEVEL
  90. #define    VISIBLE_LEVEL    80
  91. #endif
  92.  
  93. #define    Mean_Scale    2.5
  94. #define    Rate_Width    3
  95. #define    BasicSize    16
  96.  
  97. typedef    struct    {
  98.     int    lower_p, upper_p,
  99.         peak_p, peak_v,
  100.         weight;
  101.     } W_Spectrum;
  102.  
  103. typedef    struct    {
  104.     int    nsp, objects,
  105.         hillfoot_w, min_weight, height_rate,
  106.         l_u, lp, rp, weight,
  107.         slope[Rate_Width<<1], sample_w;
  108.         /*    lower, upper, width, p;    */
  109.     W_Spectrum    *s_array;
  110.     } Info;
  111.  
  112. typedef    struct    {    /* u is used different from mainpeak    */
  113.     int    l, m, u,/* any value lower than u will not be changed    */
  114.         overf;
  115.     } Threshold;
  116.  
  117. #ifdef    _DEBUG_
  118. extern    int    debug;
  119. #endif
  120. Info    info;
  121. Threshold    threshold={0};
  122. bool    adj, corner, smooth, eEn,
  123.     Msg, Bflag, Zf=16, Test, Top, idpn, errff;
  124. MType    align, frmp, fsize;
  125. int    *hist, maxhisv=256, slope_sets=Rate_Width, Valley, VP,
  126.     Ay, Ax, Ah, Aw=BasicSize, frm, bfrm;
  127. float    *mean, sum, top, enh, ovff=1;
  128. #define    Mean    *mean
  129.  
  130. #define    GValue()    arget(argc, argv, &i, &frmp)
  131.  
  132.  
  133. main(argc, argv)
  134. int    argc;
  135. char **    argv;
  136. {
  137. int    i;
  138.  
  139. format_init(&uimg, IMAGE_INIT_TYPE, HIPS, -1, *argv, "S20-1");
  140.  
  141. for (i=1; i<argc; i++)
  142.     if (*argv[i] == '-')    {
  143.     frmp=1;
  144.     switch(argv[i][frmp++])
  145.     {
  146.     case 'a':    adj++;    break;
  147.     case 'A':    Ax = GValue();    Ay = GValue();
  148.             Aw = GValue();    Ah = GValue();
  149.         if (!Aw)    Aw = BasicSize;
  150.         if (!Ah)    Ah = Aw;    break;
  151.     case 'B':    Bflag++;
  152.     case 'b':    align = GValue();
  153.         if (!align)    align--;    break;
  154.     case 'c':    top = .475;    break;
  155.     case 'C':
  156.         if (GValue())
  157.             sscanf(argv[i], "%d %d", &corner, &Aw);
  158.         else    corner = Aw = BasicSize;
  159.         break;
  160. #ifdef    _DEBUG_
  161.     case 'D':    debug++;    break;
  162. #endif
  163.     case 'e':    enh = GValue();
  164.             if (!enh)    enh=Mean_Scale;    break;
  165.     case 'h':    info.height_rate = 100 - GValue();
  166.             if (info.height_rate<5 || info.height_rate>95)
  167.                 info.height_rate = 50;    break;
  168.     case 'E':    eEn = GValue();    break;
  169.     case 'F':    ovff = GValue();    break;
  170.     case 'G':    errff++;
  171.     case 'M':    Msg = GValue();
  172.             if (!Msg)    Msg--;    break;
  173.     case 'R':    slope_sets = GValue();    break;
  174.     case 'I':    idpn++;
  175.     case 'S':    smooth = GValue();
  176.             if (smooth<1 || smooth>15) smooth=2;    break;
  177.     case 's':    info.sample_w = GValue();    break;
  178.     case 'T':    Test++;    break;
  179.     case 't':    Top++;    break;
  180.     case 'l':    threshold.l = GValue();    break;
  181.     case 'm':    threshold.m = GValue();    break;
  182.     case 'n':    frm = GValue();
  183.         if (isdigit(argv[i][frmp]))    i++,    bfrm = GValue();
  184.             break;
  185.     case 'f':    threshold.u = GValue();
  186.             if (!threshold.u)    threshold.u=32;    break;
  187.     case 'v':    Valley = GValue() - 1;    break;
  188.     case 'V':    VP = GValue() - 1;    break;
  189.     case 'z':    Zf = GValue() - 1;    break;
  190.     case 'W':    info.hillfoot_w = GValue();    break;
  191.     case 'w':    info.min_weight = GValue();    break;
  192.     case 'o':if (avset(argc, argv, &i, &frmp, 1) &&
  193.             freopen(argv[i]+frmp, "wb", stdout))    break;
  194.         message("output file -- %s", argv[i]);
  195.     default:
  196. errout:        usage_n_options(usage, i, argv[i]);
  197.     }
  198.     }
  199.     else{    if ((in_fp=freopen(argv[i], "rb", stdin))==NULL)
  200.             syserr("input file -- %s", argv[i]);
  201.         strcpy(serbuf, argv[i]);
  202.     }
  203.  
  204. io_test(fileno(in_fp), goto    errout);
  205.  
  206. (*uimg.header_handle)(HEADER_READ, &uimg, 0, 0);
  207. uimg.o_form = uimg.in_form;
  208. uimg.pxl_out = uimg.pxl_in;
  209. if (frm < 1)    frm = uimg.frames;
  210. if (bfrm < 0)    bfrm = 0;
  211. fsize = (MType)cln * row;
  212.  
  213. if (frm<2 && !Test)    syserr("Not a 3D image");
  214. if (errff) {
  215.     if (!strlen(serbuf))
  216.         strcpy(serbuf, "mean_tmp");
  217.     for (i=1; i<argc; i++)
  218.         if (*argv[i]=='-')    strcat(serbuf, argv[i]);
  219.     strcat(serbuf, ".rpt");
  220.     if (!freopen(serbuf, "w", stderr))
  221.         syserr("can't reopen %s for stderr", serbuf);
  222. }
  223.  
  224. if (uimg.in_form==IFMT_SHORT)
  225.     maxhisv <<= 8;
  226. else if(uimg.in_form)    /*    other non-byte    */
  227.     maxhisv <<= 10;    /*    1 MByte memory    */
  228.  
  229. if (!info.sample_w)    info.sample_w = SampleWidth;
  230. if (!info.hillfoot_w)    info.hillfoot_w = slope_sets*SampleWidth*Mean_Scale;
  231. if (!threshold.m)
  232.     threshold.m = fsize / ((row + cln)*pxl_bytes);
  233. while (threshold.m <= threshold.l)
  234.     threshold.m += threshold.l/pxl_bytes;
  235. if (enh){
  236.     threshold.m /= enh;
  237.     if (smooth>>2)    smooth = 2;
  238. }
  239. if (!threshold.l)    threshold.l = threshold.m>>2;
  240. if (Valley<0 || VP<0)    Valley = MAX(info.sample_w * slope_sets, 12);
  241.  
  242. ibuf = nzalloc(fsize*uimg.frames, pxl_bytes, "ibuf");
  243. obuf = nzalloc(fsize, pxl_bytes, "obuf");
  244. mean = (float*)zalloc(frm+1L, (MType)sizeof(*mean), "mean");
  245. hist = (int*)zalloc(frm+1L, (MType)maxhisv*sizeof(*hist), "hist");
  246.  
  247. (*uimg.std_swif)(FI_LOAD_FILE, &uimg, uimg.load_all=uimg.frames, No);
  248. #ifdef    _DEBUG_
  249. message("bg=%d, crnr=%d, csz=%d, ceil=%f\n", align, corner, Aw, top);
  250. #endif
  251.  
  252. {
  253. char    mesgbuf[1024];
  254. sprintf(mesgbuf, "L_threshold=%d, M_threshold=%d, fg=%d, max_v=%d, Hr=%d%%,\
  255.  MinWt=%d\nHill_W=%d, Step_W=%d, Smooth=%d, Enhanced=%.2f\n",
  256.     threshold.l, threshold.m, threshold.u, maxhisv, info.height_rate,
  257.     info.min_weight, info.hillfoot_w, info.sample_w, smooth, enh);
  258. msg("%s", mesgbuf);
  259. (*uimg.header_handle)(ADD_DESC, &uimg, mesgbuf);
  260. }
  261. (*uimg.header_handle)(HEADER_WRITE, &uimg, argc, argv, True);
  262.  
  263. if (bfrm)    fwrite(ibuf, bfrm * pxl_bytes, fsize, out_fp);
  264. i = bfrm * pxl_bytes * fsize;
  265. #define    i_buf    (char*)ibuf + i
  266. switch(uimg.in_form)
  267. {    case IFMT_BYTE:
  268.         b_get_mean_hist(i_buf);
  269.         b_align(i_buf, obuf);
  270.         break;
  271.     case IFMT_SHORT:
  272.         s_get_mean_hist(i_buf);
  273.         s_align(i_buf, obuf);
  274.         break;
  275.     case IFMT_LONG:
  276.         i_get_mean_hist(i_buf);
  277.         i_align(i_buf, obuf);
  278.         break;
  279.     case IFMT_FLOAT:
  280.         f_get_mean_hist(i_buf);
  281.         f_align(i_buf, obuf);
  282.         break;
  283.     default:syserr("unaccept format");
  284. }
  285. if (bfrm+frm < uimg.frames)
  286.     i = (bfrm+frm) * pxl_bytes,
  287.     fwrite((char*)ibuf + i * fsize, uimg.frames-i, fsize, out_fp);
  288. message("%s to %f\n", *argv, Mean);
  289. #ifdef    _DEBUG_
  290. if (debug)
  291.     for (i=bfrm, frm+=i; i++ < frm;){
  292.         message("mean(%2d) = %-f    ", i, mean[i]);
  293.         if (!(i % 3))    mesg("\n");
  294.     }
  295. #endif
  296. exit(Mean);
  297. }
  298.  
  299.  
  300. /*==============================================*
  301. *    calculate entire & each frame mean    *
  302. *==============================================*/
  303. b_get_mean_hist(buf)
  304. register byte*    buf;
  305. {
  306. register int    i, *hp0, *hp;
  307. hp = hp0 = hist;
  308. if (corner || Ax || Ay)    {
  309.     if (Ax || Ay)
  310.         message("Area = %d %d (%d x %d)\n", Ax, Ay, Aw, Ah);
  311.     for (i=sum=0; i<frm; i++, hp+=maxhisv, buf += fsize)
  312.         if (corner)
  313.         sum += (mean[i+1] = corner_bg_sample(buf, hp, corner, Aw));
  314.         else {
  315.         int    vp=0, min, max;
  316.         register int    r, c;
  317.         register byte*    bp = buf + Ay*cln + Ax;
  318.             for (r=Ah; r--; bp+=cln)
  319.                 for (c=Aw; c--;)
  320.                 hp[bp[c]]++;
  321.             find_mostvp(hp, &vp, &min, &max, maxhisv);
  322.             sum += (mean[i+1] = vp);
  323.         }
  324.     Mean = sum / i;
  325.     return;
  326. }
  327. for (frmp=sum=0; frmp < frm; frmp++)
  328.    {    register int    itmp;
  329.     hp += maxhisv;
  330.     for (i=itmp=0; i<fsize; i++){
  331.         hp0[*buf]++;    hp[*buf]++;
  332.         itmp += *buf++;
  333.     }
  334.     sum += (mean[frmp+1] = (float)itmp / fsize);
  335.    }
  336. Mean = sum / frmp;
  337. if (align)    get_sample(THRESHOLD);
  338. }
  339.  
  340. s_get_mean_hist(buf)
  341. register short*    buf;
  342. {
  343. register int    i, itmp, *hp0, *hp;
  344. hp = hp0 = hist;
  345. for (frmp=sum=0; frmp < frm; frmp++)
  346.    {    hp += maxhisv;
  347.     for (i=itmp=0; i<fsize; i++){
  348.         hp0[*buf]++;    hp[*buf]++;
  349.         itmp += *buf++;
  350.     }
  351.     sum += (mean[frmp+1] = (float)itmp / fsize);
  352.    }
  353. Mean = sum / frmp;
  354. if (align)    get_sample(THRESHOLD);
  355. }
  356.  
  357. i_get_mean_hist(buf)
  358. register int*    buf;
  359. {
  360. register int    i, itmp, *hp0, *hp;
  361. hp = hp0 = hist;
  362. for (frmp=sum=0; frmp < frm; frmp++)
  363.    {    hp += maxhisv;
  364.     for (i=itmp=0; i<fsize; i++){
  365.         hp0[*buf]++;    hp[*buf]++;
  366.         itmp += *buf++;
  367.     }
  368.     sum += (mean[frmp+1] = (float)itmp / fsize);
  369.    }
  370. Mean = sum / frmp;
  371. if (align)    get_sample(THRESHOLD);
  372. }
  373.  
  374. f_get_mean_hist(buf)
  375. register float*    buf;
  376. {
  377. register int    i;
  378. register float    ftmp;
  379. if (align)
  380.    msg("warning: No others applied for FP. Using scale to integer before %d\n",
  381.     Progname);
  382. for (frmp=sum=0; frmp < frm; frmp++)
  383.    {    for (i=ftmp=0; i<fsize; i++)    ftmp += *buf++;
  384.     sum += (mean[frmp+1] = ftmp / fsize);
  385.    }
  386. Mean = sum / frmp;
  387. }
  388.  
  389. /*==============================================*
  390. *    make all frames have same mean value    *
  391. *    or peak position or center of gravity    *
  392. *==============================================*/
  393. b_align(ibp, bbuf)
  394. register byte    *ibp;
  395. byte    *bbuf;
  396. {
  397. register int    i, itmp, diffmean;
  398. register byte    *obp;
  399. for (frmp=0; frmp < frm; frmp++)
  400.     {    obp = bbuf;
  401.     diffmean = Mean - mean[frmp+1] + top;
  402.     if (!diffmean){
  403.         obp = ibp;    ibp += fsize;
  404.     if(Msg)    message("no change for frame %-3d[%d]\n", frmp+1, diffmean);
  405.     }
  406.     else {
  407.     if(Msg)    message("mean difference(%d) = %d\n", frmp+1, diffmean);
  408.         for (i=0; i<fsize; i++)
  409.         {
  410.             if (*ibp < threshold.u){
  411.             *obp++ = *ibp++;
  412.             continue;
  413.             }
  414.             else{
  415.             itmp = diffmean + *ibp++;
  416.             if (itmp > 255)    *obp = 255/ovff;
  417.             else if (itmp < 0)    *obp = 0;
  418.             else    *obp = itmp;
  419.             obp++;
  420.             }
  421.         }
  422.         obp =bbuf;
  423.     }
  424.     i = fwrite(obp, pxl_bytes, fsize, out_fp);
  425.     if (i != fsize)    syserr("write mean data %d", i);
  426.     }
  427. }
  428.  
  429. s_align(ibp, sbuf)
  430. register unsigned short    *ibp;
  431. unsigned short    *sbuf;
  432. {
  433. register MType    i, itmp, diffmean;
  434. register unsigned short    *obp;
  435. for (frmp=0; frmp < frm; frmp++)
  436.     {    obp = sbuf;
  437.     diffmean = Mean - mean[frmp+1] + top;
  438.     if (!diffmean){
  439.         obp = ibp;    ibp += fsize;
  440.     if(Msg)    message("no change for frame %d\n", frmp+1);
  441.     }
  442.     else {    for (i=0; i<fsize; i++)
  443.         {
  444.             if (*ibp < threshold.u){
  445.             *obp++ = *ibp++;
  446.             continue;
  447.             }
  448.             else{
  449.             itmp = diffmean + *ibp++;
  450.             if (itmp > 65535)    *obp = 65535/ovff;
  451.             else if (itmp < 0)    *obp = 0;
  452.             else    *obp = itmp;
  453.             obp++;
  454.             }
  455.         }
  456.         obp = sbuf;
  457.     }
  458.     if (fwrite(obp, pxl_bytes, fsize, out_fp) != fsize)
  459.         syserr("write short");
  460.     }
  461. }
  462.  
  463. i_align(ibp, buf)
  464. register int    *ibp;
  465. int    *buf;
  466. {
  467. register int    i, diffmean, *obp;
  468. for (frmp=0; frmp < frm; frmp++)
  469.     {    diffmean = Mean - mean[frmp+1] + top;
  470.     obp = buf;
  471.     if (!diffmean){
  472.         obp = ibp;    ibp += fsize;
  473.     if(Msg)    message("no change for frame %d\n", frmp+1);
  474.     }
  475.     else{    for (i=0; i<fsize; i++)
  476.             *obp++ = diffmean + *ibp++;
  477.         obp = buf;
  478.     }
  479.     if (fwrite(buf, pxl_bytes, fsize, out_fp) != fsize)
  480.         syserr("write integer");
  481.     }
  482. }
  483.  
  484. f_align(ibp, fbuf)
  485. register float    *ibp;
  486. float    *fbuf;
  487. {
  488. register int    i;
  489. register float    diffmean, *obp;
  490. for (frmp=0; frmp < frm; frmp++)
  491.     {    obp = fbuf;
  492.     diffmean = Mean - mean[frmp+1] + top;
  493.     if (!diffmean){
  494.         obp = ibp;    ibp += fsize;
  495.     if(Msg)    message("no change for frame %d\n", frmp+1);
  496.     }
  497.     else{    for (i=0; i<fsize; i++)
  498.             *obp++ = diffmean + *ibp++;
  499.         obp = fbuf;
  500.     }
  501.     if (fwrite(obp, pxl_bytes, fsize, out_fp) != fsize)
  502.         syserr("write float");
  503.     }
  504. }
  505.  
  506. /*======================================*
  507. %    find the most frequency value h    %
  508. %    and its position (*mostvp), and    %
  509. %    effective hill feet (min, max)    %
  510. *======================================*/
  511. find_mostvp(hp, mostvp, min, max, n)
  512. register int    *hp, *mostvp;
  513. int    *min, *max, n;
  514. {
  515. register int    h, i = *mostvp;
  516.  
  517. while (!hp[i])    i++;
  518. h = hp[i],    *min = i;
  519. do if (hp[i] > h)
  520.     h = hp[*mostvp = i];
  521. while(++i < n);
  522.  
  523. while(!hp[--i]);
  524. *max = i;
  525. #ifdef    _DEBUG_
  526. if(Msg)    message("min=%d, max=%d, vp=%d, v=%d\n", *min, *max, *mostvp, h);
  527. #endif
  528. return    h;
  529. }
  530.  
  531. get_sample(T)
  532. {
  533. int    *hp=hist, f;
  534.     if (Bflag){
  535.     int    t;
  536.  
  537.     info.s_array = zalloc((MType)sizeof(*info.s_array), (MType)MOSTHILLS);
  538.     for (f=0; f<=frm; f++){
  539.     register int    i;
  540.         if (f){
  541.             if (Valley)
  542.                 for (i=0; i<Valley; i++)
  543.                 hp[VP+i] = 0;
  544.             hills(hp, maxhisv, (MType)MOSTHILLS, &info, &threshold, smooth, f);
  545.         }
  546.         else    {
  547.         int    s = (smooth!=0) + (frm<8);
  548. #if    defined    sparc | defined SUN_WORKStation
  549.         Threshold    TH;
  550.             TH.l = threshold.l * frm;
  551.             TH.m = threshold.m * frm;
  552. #else
  553.         Threshold TH = threshold;
  554.             TH.l *= frm;
  555.             TH.m *= frm;
  556. #endif
  557.             hills(hp, maxhisv, (MType)MOSTHILLS, &info, &TH, s, f);
  558.             if(!VP)    VP = info.s_array[0].upper_p - (Valley<<1)/3;
  559.             if(Valley) message("Dig valley(%d) at %d\n",Valley,VP);
  560.         }
  561.         if (info.nsp){
  562.         register int    Width=0, Weight=0, weight=0;
  563.         /*    find a hill which satisfies hill foot width    */
  564.             for (i=info.nsp; i-- &&
  565.                 (Width<info.hillfoot_w || Weight<info.min_weight);){
  566.                 Weight = info.s_array[i].weight;
  567.                 info.rp = info.s_array[i].upper_p;
  568.                 info.lp = info.s_array[i].lower_p;
  569.                 Width = info.rp - info.lp;
  570.             }
  571.             t = i+1;
  572.             if (info.height_rate){
  573.                 i=info.s_array[i].peak_v*info.height_rate/100;
  574.                 while (hp[info.lp]<i)
  575.                     info.lp++;
  576.                 while (hp[info.rp]<i)
  577.                     info.rp--;
  578.                 Weight=0;    /* recalc total weight    */
  579.                 for (i=info.lp; i<=info.rp; i++)
  580.                     Weight += hp[i];
  581.             }
  582.         /*    calc center of weight    */
  583.             Weight >>= 1;    weight = 0;
  584.             for (i=info.lp; weight < Weight; i++)
  585.                 weight += hp[i];
  586.             mean[f] = --i;
  587.         if(Msg)    message("WtCnt(%d)=%d, hw=%d, lp=%d, wp=%d, up=%d, hp=%u\n\n",
  588.             t, Weight, weight, info.lp, i, info.rp, hp);
  589.         }
  590.         else    msg("no hill found in frame %d\n", f);
  591.         hp += maxhisv;
  592.     }
  593.     if (align > VISIBLE_LEVEL)
  594.         Mean = align;
  595.     }else{
  596.     int    min, max, mostvp=T, mp;
  597.  
  598.     /*    get entire background value    */
  599.     find_mostvp(hp, &mostvp, &min, &max, maxhisv);
  600.  
  601.     /*    get bg value for each frame    */
  602.     for(f=0; f<frm; f++)
  603.        {    hp += maxhisv;
  604.         mp = T;
  605.         find_mostvp(hp, &mp, &min, &max, maxhisv);
  606.         mean[f+1] = mp;
  607.        }
  608.  
  609.     if (align > 0)    mostvp = align;
  610.     align=Mean;    Mean = mostvp;
  611.     if (mostvp > 192 || abs(mostvp-align) > 32){
  612.     msg("warning: may not properly scale to background\nentire BG = %d",
  613.         mostvp);
  614.     msg("    mean=%d : try using -C \"margin sample_size\" or -w\n",align);
  615.         for (f=0; f<frm;){
  616.             message("f%-4dBG = %-8.3f", f+1, mean[f+1]);
  617.             if (!(++f % 4))    mesg("\n");
  618.         }
  619.         return    0;
  620.     }
  621. #    ifdef    _DEBUG_
  622.     if (Msg){
  623.         msg("entire most value is %d\n", mostvp);
  624.         for (f=0; f<frm;){
  625.             msg("frame %-3d bg=%-04.0f", f+1, mean[f+1]);
  626.             if (!(++f % 4))    mesg("\n");
  627.         }
  628.         mesg("\n");
  629.     }
  630. #    endif
  631.     return    min;
  632.     }
  633. }
  634.  
  635. /*    collect given area pixel and compute mean as background value    */
  636. corner_bg_sample(buf, hp, margin, sam_sz)
  637. register byte*    buf;
  638. register int    *hp;
  639. int    margin, sam_sz;
  640. {
  641. int    i, j;
  642.  
  643.     buf += (cln+1) * margin;
  644.     for (i=0; i<sam_sz; i++) {
  645.     for (j=0; j<sam_sz; j++)
  646.         hp[*buf++]++;
  647.     buf += cln - ((margin + sam_sz) << 1);
  648.     for (j=0; j<sam_sz; j++)
  649.         hp[*buf++]++;
  650.     buf += margin << 1;
  651.     }
  652.     buf += (row - ((margin + sam_sz) << 1)) * cln;
  653.     for (i=0; i<sam_sz; i++) {
  654.     for (j=0; j<sam_sz; j++)
  655.         hp[*buf++]++;
  656.     buf += cln - margin - sam_sz;
  657.     for (j=0; j<sam_sz; j++)
  658.         hp[*buf++]++;
  659.     buf += margin << 1;
  660.     }
  661.     margin=0;
  662.     find_mostvp(hp, &margin, &i, &j, maxhisv);
  663. return    margin;
  664. }
  665.  
  666. /*======================================================*
  667. *    find group of hills & their weight in a frame    *
  668. *    peak_v < pixel_number/((r+c) * bytes_per_pixel)    *
  669. *    valley => 5 is suggested (range 1 - ...)    *
  670. *======================================================*/
  671. get_slope(ra_p)
  672. int    *ra_p;
  673. {
  674. register int    i, sum;
  675. for (i=sum=0; i<slope_sets; i++, ra_p++)
  676.     sum += Sign(*ra_p);
  677. return    sum;
  678. }
  679.  
  680. hills(hbuf, bwidth, sets, ip, trshold, smooth, f)
  681. int    *hbuf, bwidth;    /*    histo buffer width    */
  682. MType    sets;    /*    total spectrums required    */
  683. Info*    ip;
  684. Threshold    *trshold;
  685. {
  686. register int    *hp=hbuf;
  687. int    leftdir, lslope, predir=1, rightdir=0, rslope=0, ap=0, rp=slope_sets,
  688.     *minp0, *minp1, *maxp, samples=0,
  689.     main_threshold = trshold->m, lower_t = trshold->l;
  690. unsigned    min=0, max=0;
  691. bool    peakflag=0;
  692.  
  693.     if (f+eEn < 0 || f+eEn > frm){
  694.     main_threshold /= 2.5;
  695.     lower_t >>= 1;
  696.     }
  697.     if (Zf>0){
  698.     register int    i;
  699.     for (i=0; i<Zf; i++)
  700.     hp[i] = 0;    /* reset zero's value    */
  701.     }
  702.     else    msg("zeros = %d\n", *hp);
  703.     if (!Top)
  704.     for(lslope=3; lslope; lslope--)
  705.         hp[bwidth - lslope] = 0;
  706.     if (smooth){    /*    SMOOTH HISTOGRAM    */
  707.     register int    i, j, s, n = (smooth << 1) + 1;
  708.     if (idpn)
  709.         for (i=smooth; i<bwidth-smooth; i++){
  710.         for (s=0, j = -smooth; j<=smooth; j++)
  711.             s += hp[i+j];
  712.         hp[i] = s / n;
  713.         }
  714.     else    {
  715.     int *tbuf = (int*)zalloc((MType)bwidth, (MType)sizeof(*tbuf), "tbuf");
  716.     register int*    thp = tbuf + smooth;
  717.         for (i=smooth; i<bwidth-smooth; i++){
  718.         for (s=0, j = -smooth; j<=smooth; j++)
  719.             s += hp[i+j];
  720.         *thp++ = s / n;
  721.         }
  722.         memcpy(hp, tbuf, bwidth*sizeof(*hp));
  723.         free(tbuf);
  724.     }
  725.     }
  726.     while (!*hp)    hp++;    /*    pass non value zone    */
  727.     minp1 = hp-1;
  728.     {    register int    i;    /* set beingning slope    */
  729.     for (i=0; i<slope_sets; i++){
  730.         ip->slope[i] = *(hp+ip->sample_w) - *hp;
  731.         if (*hp > max)    /* get approximately    */
  732.         {    max = *hp;    maxp = hp;    }
  733.         hp += ip->sample_w;
  734.     }
  735.     }
  736.     rp = slope_sets;
  737.     leftdir = get_slope(ip->slope);
  738.     if (leftdir<0)    goto    peakset;/* if down hill, set having a peak    */
  739.     for (; hp<bwidth+hbuf && ap < sets; hp++){
  740.     if (min>*hp && *hp>lower_t)
  741.     {    min = *hp;    minp0 = hp;
  742.         if (peakflag)    minp1 = minp0;
  743.     }
  744.     else if (*hp<lower_t && !rslope)    minp1 = hp;
  745.     else if (*hp>max && rslope>=0)
  746.     {    max = *hp;    maxp = hp;    }
  747.     if (++samples < ip->sample_w)
  748.         continue;    /* not sample enough for slope    */
  749.     ip->slope[rp++] = *hp - *(hp-samples);    /* interpolate slope    */
  750.     samples=0;
  751.     if (rp < (slope_sets << 1))
  752.         continue;    /* not sample enough for a group slope    */
  753.     rp = slope_sets;
  754.     memcpy(ip->slope, ip->slope+rp, rp*sizeof(*ip->slope));
  755.     rightdir = get_slope(ip->slope+rp);    /* calc current slope    */
  756.     lslope = Sign(leftdir);    rslope = Sign(rightdir);
  757.     if (lslope == rslope)    continue;
  758.     if (rslope)
  759.     if (lslope || !lslope && Sign(predir)!=rslope){
  760.         if (rightdir > leftdir){    /* valley    */
  761. ahill:            ip->s_array[ap].upper_p = minp0 - hbuf;
  762.             max=0;    /*    change valley    */
  763.            if (peakflag){
  764.             ip->s_array[ap].weight = 0;
  765.             for (maxp=ip->s_array[ap].lower_p+hbuf;maxp<=minp0;)
  766.                 ip->s_array[ap].weight += *maxp++;
  767.             ap++;
  768.             peakflag=0;
  769.            }
  770.         }
  771.         else if (max > main_threshold){    /*    reach top    */
  772. peakset:        peakflag++;
  773.             ip->s_array[ap].peak_v = max;
  774.             ip->s_array[ap].peak_p = maxp - hbuf;
  775.             ip->s_array[ap].lower_p = minp1 - hbuf;
  776.             min = *hp;    minp0 = minp1 = hp;
  777.         }
  778.     }    /* end if lslope    */
  779.     else    if (rslope > 0)
  780.             max = 0;
  781.         else    min = *hp;
  782.     if (leftdir)    predir = leftdir;/* common exit after rslope    */
  783.     leftdir = rightdir;
  784.     }
  785.     if (peakflag)    goto    ahill;
  786.     hills_dump(ip->s_array, ap, Msg);
  787.     if(Msg){
  788.     if (hp < bwidth+hbuf)
  789.         msg("There are more than %d", ap);
  790.     else    msg("Total %d Distinguishable", ap);
  791.     msg(" Peaks in (%d)\n", f);
  792.     }
  793.     ip->nsp = ap;
  794. }
  795.  
  796. hills_dump(array, nsp, nMsg)
  797. W_Spectrum*    array;
  798. unsigned    nsp, nMsg;
  799. {
  800. register unsigned    i, j=MIN(nsp, nMsg);
  801. for (i=0; i<j; i++, array++)
  802.     message("PEAK%02d: lp->%-4d: rp->%-4d|=| pp->%-4d: pv->%-4d, w=%d\n", i+1,
  803.         array->lower_p, array->upper_p, array->peak_p, array->peak_v,
  804.         array->weight);
  805. }
  806.